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1. Introduction 

In time-dependent Hamiltonian systems the energy is generally not conserved during 
the time evolution. However, at any given moment of time it has a precise meaning. 
Namely, it is the value of the Hamilton function calculated for the particular solution of 
the equations of motion at that moment of time. Thus the system's energy depends both 
on initial conditions and on time (in the case of time-independent Hamiltonian systems 
it is specified by the initial conditions only and is conserved). In this paper we restrict 
ourselves to the consideration of the one-dimensional time-dependent harmonic oscillator 
subjected to an external forcing. This model arises in a variety of physical problems 
starting from parametric resonance jl] and ending with the geometrical ion trap quantum 
computations |2 E|- The one- dimensional harmonic oscillator was considered in a 
number of works mainly in the context of preservation of adiabatic invariants. After 
the paper by Einstein |1] , Kulsrud was the first who showed, using WKB-type method, 
that for a finite time span T the action variable / is preserved to all orders for the 
harmonic oscillator, if all derivatives of the oscillator frequency uj{t) vanish at the 
beginning and at the end of the time interval, whilst if there is a discontinuity in 
one of the derivatives he estimated the corresponding change of the action A/ [Hj. 
The same result was independently obtained for a charged particle in a slowly varying 
magnetic field for infinite time domain in [H]. Explicit general expressions were recently 
derived in |7j. Kruskal, as reported in |HI, and Lenard studied more general systems, 
whilst Gardner used the classical Hamiltonian perturbation theory Courant and 
Snyder have studied the stability of the synchrotron and analysed adiabatic invariant 
employing the transition matrix JU]- The interest then shifted to the infinite time 
domain. Littlewood showed for the time-dependent harmonic oscillator that if uj{t) is an 
analytic function, then the action I is preserved to all orders of the adiabatic parameter 
e = 1/T. Kruskal developed the asymptotic theory of Hamiltonian and other systems 
with all solutions nearly periodic [T^. Lewis, using Kruskal's method, discovered a 
connection between the action variable / of the one-dimensional harmonic oscillator 
and another nonlinear differential equation J2] • Symon used Lewis'es result to calculate 
the (canonical) ensemble average of I and its variance ^31- Finally, Knorr and Pfirsch 
proved that A/ ~ exp(— c/e) |14J. Meyer relaxed some conditions and calculated the 
constant c P^ll6j . Exponential preservation of the adiabatic invariant / for an analytic 
uj{t) on (— oo, +oo) with constant limits at t — > ±oo is thus well established Energy 
evolution in a one- dimensional time-dependent harmonic oscillator without forcing was 
studied by Robnik and Romanovski in [ZlQHl- The exact energy distribution function 
in this case was also derived [19ij . 

In this paper we consider one-dimensional classical time-dependent harmonic 
oscillator driven by an arbitrary time-dependent force. The plan is as follows. 
In section |21 exact analytical solution of the corresponding equations of motion is 
constructed in the framework of the previously developed technique based on WKB 
approach [20]. Energy evolution for the ensemble of uniformly distributed (w.r.t. the 
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canonical angle) initial conditions on the initial invariant torus is studied in the section|21 
In the same section the exact expressions for the energy moments of arbitrary order 
taken at an arbitrary moment of time are analytically derived and the corresponding 
characteristic function is analytically constructed in the form of infinite series. The limit 
of small initial ensemble's energy is considered in section |31 A simple formula coinciding 
with the one obtained in JH] for the energy distribution function is analytically derived 
in this limit. In section El the energy distribution characteristic function is numericaly 
evaluated for certain values of the system's parameters. Also, the energy distribution 
function is numerically calculated in some particular cases. Two appendices are added 
for the sake of the selfconsistency of the paper. 

We also would like to note that for the quantized time-dependent harmonic 
oscillator with an external forcing the time evolution of the displacement parameter 
can be described in terms of the classical effective Hamiltonian having the same form as 
the one studied in the rest of this paper. The particular implementation of the quantized 
time-dependent harmonic oscillator with the external forcing is the key model exploited 
in geometrical ion trap quantum computations 



2. Equations of motion and their solution using WKB approach 

We consider the time-dependent harmonic oscillator with the external time-dependent 
forcing described by the Hamiltonian 

H = ^p' + ^u;'{t)q' + ^qf{t), (1) 

where q and p are canonically conjugate generalized coordinate and momentum, uj{t) is 
the time-dependent oscillator frequency, f(t) is the time-dependent external force and 
7 is the coupling strength. Corresponding Hamilton equations of motion are 



dH 
^p 
dH 



q = ^=p, 
op 



p=-^ = -uj\t)q-^f{t). (2) 

Here the dot denotes the time derivative. This system of the first order differential 
equations is equivalent to the single second order differential equation: 

q + uj'it)q = -^f{t). (3) 

This is a nonhomogeneous ordinary differential equation and it is known that its general 
solution is the sum of the general solution of the corresponding homogeneous equation 
and some particular solution of the nonhomogeneous equation [21]. The exact solution 
of the homogeneous equation in terms of WKB series can be found in the Refs. fSl 120] 
and the summary is given in [Appendix A The general solution of the nonhomogeneous 



equation can be expressed in terms of Green function [21]. Namely, 
g(A) = qh{X) + fq = w+exp |^o-+(A)| + w_ exp |^(t_(A)| 



Energy evolution in time- dependent harmonic oscillator with the arbitrary ... 4 
J G{X,x)f{x)dx, Ao < A < A/, (4) 

Ao 

where the notations are explained in [Appendix A and fq is introduced for further 
convenience. The Green function is given by the expression [211 9.3-9]: 

Here g+ and g_ denote two hnearly independent solutions of the homogeneous equation: 
g+ = exp |^(T+|, g_ = exp|^a_| (6) 
and U{z) is the unit step function defined as 

The denominator in @ is just the Wronskian of the relevant equation and is constant 
for all X. Using the results of the [Appendix A the expression ^ can be reduced to 

G(A.x) = ^exp{r<Mj„„(£(M)f;(A_,,, (s) 

where 

A A 

r(A,x) = J A{z)dz, s{X,x) = J B{z)dz. (9) 

X X 

The requirement of the phase space volume conservation leads to the relation 
exp (^) _ 1 
B{x) ^B{\)B{x) 
Thus the Green function can be represented in the form: 



(10) 



G{\x) = ^ ^ sin(^^^V(A-x). (11) 



lB{X)B{x) 

The momentum evolution can be determined using the relationship p{X) = eg' (A) and 
it is given by the expression 

p{X) = w+a'_^{X) exp |-cr+(A)| + W-(j'_{X) exp |-(t_(A)| 

-67 J G",(A, x)f{x)dx = ph{X) + U, (12) 

Ao 

where Ph{X) is the solution of the equations of motion when the external force is zero 
(fit) = 0) and fp is proportional to the integral of Green function's derivative multiplied 
by the force. 
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3. Energy distribution moments and characteristic function 

In this section we analyticaUy derive the exact expressions for the energy distribution 
moments taken at arbitrary time t = ti for the ensemble of uniformly distributed (w.r.t. 
the canonical angle) initial conditions on the initial invariant torus. Using them we 
construct the characteristic function of the energy distribution in form of an infinite 
power series. First, we note that the generalized coordinate and momentum of the 
system at time ti are the linear combinations of the same quantities taken at the initial 
time moment t = tQ plus nonhomogeneous terms originating from the forcing. Namely, 

where qi = q(ti), pi = p(ti), go = <l(to), po = p(to), the coefficients a, b, c and d are 
the same as in the homogeneous case and the corresponding expressions can be found 
in the [Appendix A The quantities fq and fp are defined in (jH) and (fT^ and are given 



by the expressions 

fq = -1 G{\,x)f{x)dx, 



Ao 

Af 



U = -e^ G'^{\x)f{x)dx. (14) 



Ao 

The energy of the system Ei = H{q(ti),p{ti),ti) at the time t = ti is 

Ei = \pl + \ujlql + ^q^h, (15) 

where uJi = u!{ti) and /i = /(ti). If the initial energy of the ensemble equals 
Eq = H{q(to),p(tQ),to) then using (fT^ . assuming /o = (without losing generality) 
and performing the transformation 

^o = J-^cos0, Po = j2E'osm(f) (16) 
V ^0 

with ujq = a; (to), one obtains 

El = a cos^ <P + P sin^ + 7 sin 20 + a cos + p sin + 5. (17) 



Here 



Eq / 2 2 2 

a = ^ (c +uj,a 
(3 = Eo (d^ + i^lH' 
7 = — ^ (cd + uolab 



l2Eo 







cr = J ^ (c/p + ulafq + 7a/: 



1 ' 



p = pEo(dU + ulbfq + ^bU), 

1 f2 I ^, ,2 f2 
2 



/p + o^r/.+7/Ji. (18) 
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Using the equation (jl7|) one can calculate the average energy of the ensemble at time 

E^ = ]^{a + l3) + 5 (19) 
and the deviation of the energy from its mean value 



Ei-Ei = + 72 sin (20 + V-i) + V^' + sin (0 + ^2), (20) 

where 

K = - (a — /?) , tan-?/'! = — , tajitp2 = —. (21) 

2 7 p 

Using equation (j^n|) we obtain the exact expressions for the even and odd moments 
of the energy distribution. The main steps of a bit bulky derivation are given in 
the [Appendix B The results are 



X ( ~ M COS [(2n - 2p - 2g + 1) ^1 - (2p - 2s - 1) ^/-a] 



s 



n— 1 n— p— 1 p— 1 



x52(2n-2p-2q+i),2p-2s-i + X! X! X! — ^i;^::! — ^^-B^^ 

p=l g=0 s=0 ^ 

^2^)1 ) ( ) - 12 (" - p - - 2 (p - .) *d 

El - ^i) = V V V I ^2n-2p^2p-l 2n 1 

X ( ^^"""^^ ^ ^^"^ ^sin[2(n-p-g)^i-(2p-2s-l)^2] 

x54(„-p-.),2p-2.-i + E E E ^4^r^^""'"^'' 

p=l g=0 s=0 

^ 1^ 2n - 1 ^ 1^ 2(n - j9) - 1 ^ 1^ 2j9 ^ 

X sin [(2n - 2p - 2g - 1) V^i - 2 (p - s) ?/'2]52n-2p-2q-i,p-s- (23) 



Here 



A = + 72, 5 = ^cr2 + p2. (24) 
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For any specific numerical value of n the sums in (I22p -(|23 p can be analytically evaluated 
using one of the algebraic software packages, for example Mathematica. Below we 
provide the expressions for the few lowest moments: 

(^1 -E^f = -^AB^ sin (V^i - 2^2), 

{E,-E,y = 1{a' + 4A'B' + B') . (25) 

The first moment (^Ei — E^ = as it should be. Equations ((221)- (ESI allow us to 
construct the characteristic function of the energy distribution in the form of an infinite 
series according to the following formula: 

It is seen that the characteristic function is complex. It means that the energy 
distribution is not an even function of Ei — Ei in contrast to the case of the time- 
dependent harmonic oscillator without forcing |7l |^ 120] • The same conclusion 
follows from the equation (|^^. since the odd moments are non-zero in the general case. 



4. Small initial energy limit 

The expressions (j22j) - (j23|) can be significantly simplified if we consider the evolution 
of the initial ensemble distributed close enough to the point (g = 0,p = 0) and thus 
having small initial energy Eq. In this case the series for the characteristic function ()26|) 
can be summed up and the expression for the distribution function can be analytically 
obtained. From equations (fTH|) . ()24|) it follows that A Eq and B ~ ^Eq , where 
^ has the dimension of the square root of the energy, but it does not depend on Eq. 
Therefore when Eq is small, A B. Starting from either (j^UI) or directly from (j^^ -()23 p 
it follows that the even moments in this case are the same as for the time-dependent 
harmonic oscillator without forcing considered in the Refs. ^1 ^] and are given by 
the expression 



-V2^_ (2n-l)!! 



'E^-E,y^^ '"'^ ^"- /x^", (27) 

where /x^ ~ B^/2 is the second moment or variance. The odd moments are zero in this 
approximation. The average energy value equals Ei ^ S and has also the meaning of 
the energy acquired by the system at the time ti starting its motion at the initial point 
(go = 0,po = 0). The characteristic function in this case is given by [TU] 

00 1 / 2 2 \ 
/(■!')«E(^(-^) =^«(«V2). (28) 
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Figure 1. The plots of the real and imaginary parts of the energy distribution 
characteristic function for the certain values of the parameters. Namely, A — 0.7, 
B = 0.7, V'l = 7r/4 and ■(/'2 = 7r/4. 



where Jo{z) is the zero order Bessel function. It yields the following distribution 
function [T^ : 

P(Ei) = , ^ (29) 

7rW2/i2 - (e, - El" 



This function has two singularities at (^Ei — Eij = ±/iA/2 and it is even with respect to 
El — El). 



5. Numerical results 



The characteristic function in ()26p can be numerically evaluated for certain values of 
the parameters. The plots of its real and imaginary parts, when A = 0.7, B = 0.7, 
ipi = 7r/4 and 1^2 = 7r/4, are presented in figure^ Here the first forty energy distribution 
moments have been summed up in order to approximate the characteristic function. We 
also extended the summation up to the sixty moments (not shown) in order to convince 
ourselves that the result does not change. 

Now we numerically obtain the energy distribution function. For this purpose 
we exploit equation ()20|) . We divide the range from to 2tt in 10000 intervals, let 
have a value in each of this small intervals and find how many corresponding values of 
X = (^Ei — Ei^ fall into one of the 200 intervals covering the range from — (^A + 
to A + B. After the normalization we obtain the histogram approximating the energy 
distribution function. We perform the corresponding computations for several values of 
the parameters. Primarily, we considered the case when A = 0.7, B = 0.7, ipi = tt/A 
and '?/'2 = 7r/4 (the same values as used in the figure ^ . The result is shown in figure El 
It is seen that the distribution function has four singularities in this case instead of two 
reported for the time-dependent harmonic oscillator without external forcing [3 
or obtained in the small initial energy limit in section |3] Two additional singularities 
appeared due to the additional term in proportional to sin (0 + '?/'2)- Also the 
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Figure 2. The normalized histogram approximating the energy distribution function 
P{x), when A = 0.7, B = 0.7, Vi = 71-/4 and V'2 = 7r/4. 



0.6- 



0.4- 



-0.5 



. 5 



Figure 3. The normahzed histogram approximating the energy distribution function 
P(x), when A = 0.5, B = 0.5, ipi = tt/2 and ip2 — vr/4. The distribution fmiction is 
symmetric with respect to x. 



distribution function is not even with respect to x as expected. These conclusions are 
supported by the numerical calculations at other values of the parameters (not shown). 
However, we found that the distribution function can be a symmetric one with respect 
to X for some particular parameter values. For example, this is the case ii A = 0.5, 
B = 0.5, ipi = 7T /2 and iIj2 = 7r/4, where according to the third moment is zero. 
The result of the computation is shown in figure El Indeed, it is clear from (j^n|) that 
P{x) is an even function of x if ipi = 2ip2 in which case all odd moments ()2Hj) vanish, 
which has been confirmed using Mathematica up to and including 79th moment. We 
have also checked the behaviour of the distribution function in the small initial energy 
limit. The corresponding plot is given in figure lU The following parameters values 
were chosen: A = 0.007, B = 0.7, ipi = and ip2 = It is seen that the energy 
distribution function has only two singularities in this limit as it should be according to 
the equation and it is symmetric w.r.t. x. 
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Figure 4. The normalized histogram approximating the energy distribution function 
P{x) in the small initial energy limit. Parameters values are A = 0.007, B = 0.7, 
-01 — 7t/3 and ip2 = 7r/4. The distribution function has only two singularities. 

6. Conclusions 

We have constructed the general solution of the equations of motion for the time- 
dependent harmonic oscillator with external forcing. The approach based on WKB 
method developed in [3 UHl dl 120] was exploited for this purpose. Using it we have 
analytically derived the exact expressions for the moments of the energy distribution 
of arbitrary order calculated at arbitrary moment of time for the case of the ensemble 
of uniformly distributed (w.r.t. cannonical angle) initial conditions on the invariant 
torus. We have obtained the characteristic function of the energy distribution in the 
form of an infinite series. It appears to be a complex one, since the odd moments are 
non-zero in distinction to the case of the time-dependent harmonic oscillator without 
forcing considered in [71 d]. The energy distribution function is not even w.r.t. 
X = (^Ei — Eij in the presense of external forcing. We have considered the limit of the 
small initial ensemble's energy and analytically obtained the expressions for both the 
characteristic function and the energy distribution function in this case. They coincide 
with the expressions for the same quantities in the case of the time- dependent harmonic 
oscillator without forcing but the physical meaning and origin of the parameters is 
now different. The characteristic function is real and the energy distribution is an 
even function w.r.t. x = {Ei — Ei) and has two singularities in this limit. We have 
numerically obtained the characteristic function and the energy distribution for some 
particular values of the parameters for the general case, where the energy distribution 
function has four singularities, but in the small initial energy limit their number reduces 
to two, indeed. The results obtained can be useful in the context of geometric ion trap 
quantum computations where the particular implementation of the investigated 

model of the time-dependent harmonic oscillator with external forcing plays the crucial 
role. 
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Appendix A. Solution of the homogeneous equations of motion in the 
framework of WKB approach 

In this section we give a brief summary of the results concerning the solution of the 
homogeneous equation 0, which is 

q + uj\t)q = 0. (A.l) 

More details can be found in [71 QB]- Let us first introduce the dimensionless time 
A = et, where e = 1/T and T is the typical time scale under consideration of the 
system's dynamics. Equation (jA.lll is transformed to 

eV(A)+a;2(A)g(A) = 0. (A.2) 

The prime denotes the differentiation with respect to A. The solution of the 
equation ()A.2|1 is sought in the form 

g(A) = ?i;exp|i(T(A)|. (A.3) 
Here w is some constant with the dimension of length. After the substitution one gets: 
{a'{\)f + ea"{\) = -u\\). (A.4) 
The WKB expansion for the a is 

oo 

o-(A) = 5:eV,(A). (A.5) 

k=0 

The substitution of ()A.5j) into ()A.4|1 and the comparision of the coefficients at the same 
powers of e gives the following recursion relation: 



^0 = ^1 = -TTT' 

ZCTq 

^" = ^'k^'n-k + , > 2. (A.6) 

Equation ()A.2|) has two linearly independent solutions given by the choice of a'^^^X) = 
iLj{X) or crQ_(A) = ~iuj{\) as the starting term in the recursion ()A.6|) . Since the 
frequency u;(A) is real, it follows that all functions cr'2k_^_i are real and all functions 
are pure imaginary and the following relations hold: 

'^2k,+ — ~^2k-^ ^2k+l,+ — ^2k+l- — ^2k+l^ i.^-^) 

where A; = 0, 1, 2, . . .. Thus 

a'^{\) = A{X)+iB{\), a'_{\)=A{\)~iB{\), (A.8) 

where 



oo 



A{\) = Y: e^'^V^,^,(A), B{X) = -^Y: e''a',,^^{X) (A.9) 

fc=0 fc=0 

are both real quantities. The integration yields 

a+{X) =r{X)+ts{X), cx_{X) =r{X)-ts{X), (A.IO) 
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A A 

r(A) = j A{z)dz, s(A) = J B{z)dz. (A.ll) 

Ao Ao 

Here Aq corresponds to the initial moment of time. Thus the general solution of the 
homogeneous equation ()A.2j) is given by 



qh{\) = w+exp 1^0-+ (A) I + w_exp|^(j_(A)|. 



(A.12) 



The initial values of the coordinate and momentum of the homogeneous system 
{f{t) = 0), namely, go = Q'(^o) and po = p(to) are connected with their final values 
qi = q{ti) and pi = p{ti) at time t = ti by the linear transition map: 



where 



$ : 



Po 



Pi 



a b 
c d 



(A.13) 



(A. 14) 



and the matrix elements are given by the expressions [3 ^] : 



5(Ao)5(Ai) 
1 

sin 

fi(Ao)fi(Ai) 

1 



yl(Ao) sin 
^(Ai)\ 



B{\q) cos 



'six, 



(A(Ao)A(Ai)+S(Ao)S(Ai))sin 



's(Ai 



+ (A(Ao)5(Ai)-A(Ai)5(Ao))cos 



's(Ai 



1 



5(Ao)5(Ai) 



A{Xi) sin 



's(Ai 



B{Xi) cos 



's(Ai 



(A.15) 



Appendix B. Derivation of the exact expressions for the energy moments 

In this section we provide the main steps of the energy distribution moments derivation. 
Two cases have to be considered. Namely, the case of the even moments and the case 
of the odd ones. We start with the even moments of arbitrary order. They are 

2n 1 

~ 2^ 



Ei-E,) =—jd<P[E,-E, 



1 

= — j d0{isin(20 + Vi) + 5sin (0 + ^/'2)} 



(B.i; 
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where the equation (^01) was used, and A = 1/^2 + ^2 g^^^j _ ^^-2 _|_ p2_ Applying 
Newton's binomial formula to ()B.1|) and expanding the powers of sines fIT . 1.320] the 
even moments can be represented in the form: 

- \2n 1 ~ 



El - Ei^ = — i^" J c/0sin2" (20 + ^1) 



27r f::^ \ 2p - I 



2tt 



p=l 

2n-2p+l 
Q 

1 



22n—2p 



E(-i) 



n—p+q 



q=0 



sin ([2n - 2p - 2g + 1] [20 + ^/'i]) 



22p-2 



E(-i) 



p+s-1 / 2p - 1 



1 _ _ 

+ _ ^2n-2p^2p 



( 2p ) i) 

X cos {2[n — p — q] [20 + -^i]) 



sin ([2p- 2s- 1] [0 + V'2]) 

^ n-p-l 

J2 (-1)""P"^2 

g=0 



22{n-p) 

2(ri — p) 



2{n — p) 





1 1 

22p 


"p-1 


•1 


E( 






/•27r 


+- 






271 


Jo 



2p 

s 



y n — p 

cos {2[p~s] [0 + ^/'2]) 



2p 



/ d0sin2"(0 + ^2 

JO 



(B.2) 



The integration over the angle yields the result (see 2.511]): 

27r , ^ . . 

-x2n (2n-l)!! 



/ - \2n If 



n-pp-1 ( _^yi+q-{-s-l 



+ EEE 

p=l q=0 s=0 

' 2p-l 
s 



22n- 



El — El 

^2n-2p+1^2p-l 



2"n! 

2n 
2p- 1 



(i^" + 52n 



2n - 2p + 1 
9 



X 



cos [{2n - 2p - 2g + 1) ^1 - (2p - 2s - 1) 1P2] 

n-l n-p-l p-1 ('_-|^^n-g-s 



X^2(2n-2p-2i}+l),2p-2s-l + E E E 

p=l g=0 s=0 

2n\ f 2{n-p) \j 1^ 2p 

n-l 



_^2n-2p j^2p 



X 



2j9 



x5. 



1 



2(n-p-g),p-s ' ^ f,2n 
p=l ^ 



22"- 1 

COS [2 {n — p — q) ipi — 2 {p — s) 1^2] 
2n \ I 2(n — p) \ I 2p \ 



2p 



n — p 



p 



(B.3) 



Calculation of the odd moments is performed in the same way. Namely, 

2-K 



f - \2n-l \ ( f 

{El-El) =—jd<p[Ei-Ei 



2n-l 
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f r ~ ~ -I 2n-l 

y rf0{Asin(20 + V^i) + 5sm(0 + V^2)} • (B.4) 



27r 



Applying Newton's binomial formula and expanding the powers of the sines j22| 1.320] 
we obtain 



2p-l J J ^22("-p) 



\-l)"-^-'2 ( ] COS [2 (n - J9 - g) (20 + i;,)] 



g=0 \ (1 

p-1 



+ 



1 rf^sin {{2n - 2p - 2g - 1) (20 + ^i)] — 



27r 

X 





X:(-l)^-^2 j cos [2 {p - (0 + ^2)] + ( ) ^ (B-5) 

where we took into account that the definite integral of the odd power of the sine taken 
over the sine period is zero. The remaining integrations give the result: 

n— 1 n-p— 1 p— 1 



p=l g=0 s=0 ^ \ ^f* -L / 

X ( ^^""^"^^ ^ '^^'^ ^sin[2(n-p-g)^i-(2p-2s-l)'02] 

x54(„-p-.),2p-2.-l + E E i-2 i^"-^^-^^^^ 

p=l q=0 s=0 

^ 1^ 2n - 1 ^ 1^ 2(n - - 1 ^ 1^ 2j9 ^ 

X sin [{2n - 2p - 2g - 1) - 2 (p - s) ^/'2]52n-2p-29-i,p-s- (B.6) 

Thus we have analytically derived the exact expressions for the even and odd energy 
distribution moments of arbitrary order. 
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